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The interface between the vapor and liquid phase of quadrupolar-dipolar fluids is the seat of an 
electric interfacial potential whose influence on ion solvation and distribution is not yet fully under- 
stood. To obtain further microscopic insight into water specificity we first present extensive classical 
molecular dynamics simulations of a series of model liquids with variable molecular quadrupole mo- 
ments that interpolates between SPC/E water and a purely dipolar liquid. We then pinpoint the 
essential role played by the competing multipolar contributions to the vapor-liquid and the solute- 
liquid interface potentials in determining an important ion-specific direct electrostatic contribution 
to the ionic solvation free energy beyond the dominant polarization one. Our results allow us to 
expose shortcomings in both purely dipolar models of water near interfaces and recent mesoscopic 
statistical mechanical approaches to evaluating ionic potentials of mean force near interfaces. 



I. INTRODUCTION 

At the macroscopic interface between a liquid (/) and 
its vapor (v) phase there is a spatial inhomogeneity that 
induces a charge imbalance, producing an electric field 
and consequently a potential difference across the inter- 
face, <f>i v = 4>i — (f> v . Despite extensive molecular simula- 
tion studies at both the classical and quantum mechan- 
ical levels over the past few decades fTl-TH]. a complete 
understanding of this potential, how it depends on the 
characteristics of the fluid studied, and its role in the sol- 
vation of ions is not yet at hand |15| I16j . Such an under- 
standing has become a pressing matter, because there is 
currently much interest in constructing mesoscopic mod- 
els of electrolytes near vapor-liquid interfaces and solid 
membrane surfaces and in nanopores [15-24, 26 . It is 
also now clear that the value of the interface potential 
observed experimentally depends on the probe used: al- 
though electron diffraction and holography techniques 
may measure the full interface potential, electrochemical 
techniques involved in ion solvation seemingly do not [T2T - 

mng. 

Until now mesoscopic approaches to ion distribution 
have either completely neglected the contribution of the 
interfacial potential [El l22Tf2"4] or, as we demonstrate 
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here, severely overestimated its importance by treating 
finite size solutes as point test charges [5UJ 03] ■ For fi- 
nite size ions a second microscopic solute-liquid interface 
potential, 0; s = (f>i — <p s , exists, defined as the poten- 
tial difference between the bulk liquid (I) and the center 
of a neutral solute (s) (fig. [I]) [whose size is determined 
in Molecular Dynamics (MD) simulations by the short 
range repulsion of the Lennard- Jones potential] . The po- 
tentially important role played by this microscopic poten- 
tial in determining ion distribution near inhomogeneities 
needs to be clarified (TUJ in order to provide deeper theo- 
retical insight into both molecular simulations and exper- 
imental results. Furthermore, the coupling between the 
quadrupolar and dipolar contributions to the interface 
potential and their respective roles in governing ion dis- 
tribution need to be reconsidered. To do so we present ex- 
tensive classical molecular dynamics simulations of model 
liquids that interpolate between SPC/E water (a classi- 
cal three site partial charge model [57]) and a purely 
dipolar liquid. 

In physical terms our study can be viewed as part 
of the quest, still far from complete, for the physical 
components of the position dependent ionic Potential 
of Mean Force (PMF), $(r), near dielectric interfaces 
and surfaces arising from solvent-ion and ion-ion inter- 
actions after the solvent degrees of freedom have been 
integrated out [HI H9H24] , We focus uniquely on the 
poorly understood role played by the interfacial poten- 
tial in determining the electrostatic contribution to ion 
solvation. Although other contributions to the ionic PMF 
and solvation-free energy, such as the hydrophobic and 
dispersion ones, may play non- negligible roles and there- 
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fore be important for interpreting molecular simulations 
and understanding experiments, these contributions will 
not be considered here (as they are already fairly well un- 
derstood thanks to recent progress in this area) [TM2"5] . 

One major impediment to obtaining the physically 
identifiable mesoscopic contributions to the ionic PMF, 
<I>, arises from questions concerning the amplitude and 
sign of the vapor-liquid interface potential and the role 
it plays in determining the ionic solvation free energy. 
In order to address these questions we compare a direct 
evaluation from MD simulations of the two relevant elec- 
trostatic contributions to the ionic solvation free energy 
for SPC/E water - a direct interfacial one that does not 
account for solvent polarization due to the ionic charge 
and a polarization one that does - with simplified ap- 
proaches previously adopted in the literature (namely, a 
direct one approximating ions as point test charges and 
a simple Born-type polarization approximation, defined 
below). 



II. VAPOR-LIQUID INTERFACE POTENTIAL 
AND IONIC PMF: STATE-OF-THE-ART 



Near a planar vapor- liquid interface the local ion 
concentration can be expressed in terms of the PMF, 
as Pi{z) — pu exp[— $(z)/J;bT], where z is the 
normal coordinate and pu is the ionic concentration in 
the bulk liquid (where $ is taken to vanish). The to- 
tal ion solvation free energy can then be expressed as 
AG lon = — <&(z v ), where z v is in the vapor phase (see 
fig. [I]). In theoretical studies of both vapor- liquid wa- 
ter interfaces and membrane-liquid surfaces, it has some- 
times been hypothesized [2"0H2"5] that the bare interface 
potential enter the PMF via a simple direct electrostatic 
contribution $p Qt (z) — q[5cj)(z) — where q is the ion 
charge, 5<j>{z) = 4>(z) — <f> v is the local value of the poten- 
tial difference and <f>i v = Scp(zi) (zi is at the center of the 
liquid slab far from the interface). This approach, which 
amounts to treating a finite size ion as a point test charge 
q [2"0Tl23| . is critically examined here for SPC/E water. 

Classical molecular dynamics (MD) simulations pre- 
dict potentials on the order of —0.5 V for both vapor- 
liquid interfaces and membrane-liquid surfaces and, if 
used in the point ion approximation, $ pot , would seem- 
ingly yield the dominant contribution (~ 20 k^T for 
monovalent ions) to the PMF over a substantial part of 
the interfacial region |21j . This approximation, however, 
predicts incorrect results for the PMF and correspond- 
ing ion density, when compared with MD simulations, 
both in the infinitely dilute limit |21j and when incor- 
porated into a modified Poisson-Boltzmann approach (to 
study higher electrolyte concentrations [30] ): neither the 
strong build-up of anions near a strongly hydrophobic un- 
charged surface ([20], Fig. 4a), nor the variations in the 
dilute limit of the PMF near a membrane surface (|21j. 
Fig. 3) predicted by this approach are in agreement with 
MD simulations. This approximation also yields a very 
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FIG. 1: Electric potential variation, 5(f>(z), across the vapor- 
liquid and the liquid-solute interfaces for the neutral LJ solute 
1° immersed in liquid water (SPC/E), at z = 50 A. The Gibbs 
dividing surface (GDS) of the macroscopic interface and the 
solute position are indicated by dashed vertical lines. 



substantial, albeit seemingly undetected, direct contribu- 
tion to the ionic free energy of solvation, 



AG' = -% ot (z v ) 



(i) 



on the order of 25 k^T. Disturbingly, the reasonable 
agreement between the experimental results for the sur- 
face tension of electrolyte solutions and certain promising 
mesoscopic theoretical approaches that neglect the in- 
terface potential completely would be severely disrupted 
if such large interfacial potentials were taken into ac- 
count [22j [23] • This situation becomes even more com- 
plicated if one considers that more "realistic" quantum 
mechanical calculations can lead to positive interface po- 
tentials of much higher amplitude (+3 eV) [TT]-[T3], but 
no signature of such a potential is seen in recent ab initio 
simulations of ion solvation 14J . 



III. IONIC FREE-ENERGY OF SOLVATION 

The electrostatic (ES) contributions to the total ion 
solvation free energy for the models studied here can be 
extracted directly from the MD simulations via 



AG'^i = AG + AG pol 



with 



and 



AG = 



= q{<Piv 



AG pol ~ <0% - cj> sv )/2 



(2) 



(3) 



(4) 



3 



where <j> sv — 4> s - <j> v = 4>i v - 4>i s and t^ 1 ™ = 4>i v - 0J° n are, 
respectively, the total vapor-liquid-solute potential vari- 
ations for neutral and charged solutes. The ion solute- 
liquid interface potential, ^>J° n = <^ on — 4>f n is the po- 
tential difference between the bulk liquid and the cen- 
ter of the charged ion (with the bare Coulomb potential 
due to the ion itself subtracted out). The first, direct, 
term, AGo = q4> sv , may be regarded as the electrostatic 
free energy of solvation of an ion placed in the poten- 
tial <p sv created around its corresponding neutral counter- 
part; the second polarization term, AG po i, obtained from 
an approximate generalized "charging" method [251 j 
arises from the response of the solvent to the ion (which 
generates an overall potential variation t/) 1 ™ much larger 
in amplitude than <f) sv ) [29] . Because <pi v is ion inde- 
pendent, the polarization contribution can be simplified 
to AGpoi — q{<f>i s — 4>\° n )/2 and therefore be obtained 
from bulk simulations (this contribution is the micro- 
scopic analog of the mesoscopic Born one presented be- 
low). 



A. Mesoscopic Born model 

Within the mesoscopic Born model, an ion is modeled 
as a point charge q sitting in its spherical cavity of ef- 
fective radius Ri in bulk water treated as a continuum 
of dielectric constant e w . The radial electric potential 
around the central ionic charge, l ~p- ma (r\ determines the 
Born approximation for the polarization contribution via 



AG^ = I limfe on (r) - Vo (r)] = 



q 

8Tre Q Ri 



1 



- 1 



(5) 

where eo is the vacuum permittivity, <po(r) is the bare 
Coulomb potential, and e w ~ 78 is the dielectric con- 
stant of bulk water at room temperature. Although this 
type of polarization contribution (~ 100 — 180 k^T) typi- 
cally dominates the ionic solvation free energy for the ions 
studied here, it is neither clear how accurate the simple 
Born approximation is (due to the neglect of potentially 
important ion-solvent correlations near the ion) , nor how 
to choose best R4. 

Furthermore, despite its dominant role in the global 
ionic free energy of solvation, the polarization contribu- 
tion to the local ionic PMF, $(z), is seemingly not the 
dominant contribution over a significant part of the inter- 
facial region, which means that the role of other contri- 
butions must be clarified [2T] . Although a Born- type po- 
larization term is commonly incorporated in mesoscopic 
approaches to the free energy of solvation (or PMF), the 
direct term, arising from the bare interfacial potential, is 
either completely neglected without justification [T!?II221 - 
[23] or strongly overestimated by incorrectly assuming 
4>i s = in AGo (eq. [3 1, which leads to the approxima- 
tion AG' = q <j>i v (eq. 11) for the direct contribution [or 
$ pot (z) in the PMF, which includes only the contribution 
of the vapor-liquid interface [501 HI] and neglects entirely 



that of the solute- liquid one] (see fig. 1). 



IV. MODELS AND METHODS 

We compare a direct evaluation from MD simulations 
of the two terms contributing to AG^g, namely AGo 
and AGpoi, for SPC/E water with the simplified approx- 
imations presented above, respectively, AG and AGp 0l . 
In order to shed further light on the interplay between 
the solvent molecular dipole and quadrupole moments 
(and thus water specificity), a series of molecular models 
having the same permanent dipole moment as SPC/E, 
but different quadrupolar ones, were first generated by 
reducing the H-O-H angle 7, while keeping fixed both 
the original partial charges on each site and the distance 
between the oxygen and the midpoint between hydro- 
gen atoms. For all but one molecular model the SPC/E 
parameters were maintained for the Lennard- Jones (LJ) 
interaction centered on the oxygen atom. For the n th 
molecule of each liquid model we define the molecular 
dipole moment p„ = Y] ■ qjVj n and quadrupole moment 

tensor (Q n ) a p = Qj®jn,aXjn,p, with x jn, a the a th 
Cartesian component of the position vector r jn of par- 
tial charge qj (j = 1,2,3) with respect to the center of 
charges within molecule n. We can then compute the 
macroscopic polarization 



P(r) = (£ Pn S(r-r n ) 



and the macroscopic quadrupole moment density, 
Qap(r) = g (^2 n (Qn) a /3S{r - r„) 



(6) 



(7) 



directly from the simulations as ensemble averages [2J. 

The local electric charge density, p(r), can be eval- 
uated directly by extracting the partial charge density 
associated with the particular molecular model; and the 
associated electric field E and potential <f> can then be 
obtained by integration of the Poisson equation in ap- 
propriate coordinates: 



-V-E = 



t]_ 



(8) 



The mean electric field along the direction normal to the 
planar vapor-liquid interface, at distance z is written in 
Cartesian coordinates as: 



E z {z) 



£0 



dz , 



(9) 



where p {z') represents the average electric charge den- 
sity (evaluated within the scope of MD simulations as the 
volume density of the sum of partial charges associated 
with a particular molecular model found in the bin cor- 
responding to the position z'). The coordinate z v is the 
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origin of integration in the vapor phase (far from the in- 
terface). Similarly, in spherical coordinates (appropriate 
for the curved solute-liquid interface), the mean electric 
field is given by 



H J e Q 



(10) 



The total charge density can also be expressed as a 
multipole expansion (30] . 



P = -V ■ P + V a Vf)Q a f} 

a,/3 



(11) 



here truncated after the quadrupolar term. The dipo- 
lar and quadrupolar electric fields and potentials can 
then be obtained from the dipole moment density P and 
quadrupole moment density Q, respectively: 

V 2 D = -V-E D = -V-P = --p D (12) 



v- 



£0 



-P Q (13) 



After computing the full interface potential from eq.[HJ we 
compare it with the sum of the dipolar and quadrupolar 
contributions computed from eqs . [12] and [13] to assess the 
accuracy of the truncated multipole expansion. 

The dipolar contribution to the total electric field can 
be obtained from p D = — V • P, the mean dipolar charge 
density created by the distribution of the macroscopic 
dipole moment density P, yielding for the z-component: 



£o 



(14) 



since P z vanishes at z v in the vapor phase, far from the 
interface region. Similarly, the dipolar component of the 
radial field in spherical coordinates, appropriate for the 
solute-liquid interface, reads: 



Pr(r) 



(15) 



where P r (r) represents the radial distribution of the den- 
sity of the dipole moment. The dipole moment density P 
obtained from the MD simulations as ensemble averages 
of molecular dipole moments is presented in detail below 



for the planar 1-v interface (Section VA| 



The determination of the mean electric field (or charge 
density) permits the calculation of 8(j> (z) — <p (z) — <f> (z v ), 
the local electric potential difference evaluated at posi- 
tion z in the vicinity of the interface: 



•(*) 



E z (z') dz' = — p (z r ) (z' - z) dz' (16) 



Similarly, the dipolar, S(j> D (z), local potential profiles can 
be obtained from the corresponding electric field, (z) 
(or charge density, p D ). 

The quadrupole moments of the models SPC, SP9- 
SP5, and S2N range from the SPC/E value down to zero 
(table [I]). Because both dipolar S2N and S2L models are 
asymmetric, a symmetric dumbbell-like model (S2D) was 
also investigated (with two LJ spheres on both ends of 
the dipole). Ions are modeled as simple point charges 
carrying an L J sphere. Simulations of the vapor- liquid in- 
terface were carried out using a modified parallel version 
of the molecular dynamics package Amber 9 [3Tj and a 
slab geometry methodology similar to the one often used 
in the literature: 1000 liquid molecules placed in a rect- 
angular unit cell of dimensions 31.04 x 31.04 x 91.04 A 3 
occupying roughly the middle one-third of the available 
space and generating two vapor-liquid interfaces [2H [S3] • 
A Lennard- Jones interaction potential is centered on the 
solvant oxygen atom, characterized by the SPC/E pa- 
rameters a — 3.1657 Aand e = 0.1553 kcal/mol. 

In the "bulk" ion solvation simulations, the system 
comprised a cubic cell of 1000 water molecules, with one 
solute immersed at its center. The ion properties are 
summarized in table [IT] Each ion is modeled by a simple 
point charge, a Lennard- Jones potential defined by a and 
e and, when polarizable models are analyzed, a polariz- 
ability. The cross parameters for the ion-water Lennard- 
Jones interaction are determined via Lorentz-Berthelot 
mixing rules. 

Periodic boundary conditions were applied in all 
three directions. The long-range charge-charge, charge- 
dipole and dipole-dipole interactions were treated by the 
particle-mesh Ewald summation method for both the 
charge and dipole moments [38 . For computational ef- 
ficiency, in the polarizable simulations, an extended La- 
grangian method was utilized to compute the induced 
dipole moments, regarded as additional dynamic vari- 
ables [15]. 

A cutoff radius of 10 A was used for the short-ranged 
non-bonded LJ interactions and for the real space compo- 
nent of the Ewald summation. The geometries of the liq- 
uid molecules were constrained by applying the SHAKE 
algorithm with a relative geometric tolerance of 10~ 4 . 
The equations of motion were integrated using the veloc- 
ity Verlet algorithm with a default time step of 1 fs [I0J . 
In order to avoid occasional drifts of the slab along the z- 
axis normal to the interface, the center-of-mass (COM) 
velocity was removed every 1000 steps. Configurations 
were saved every 100 fs in the output trajectories and 
each such frame was readjusted with respect to the z- 
axis, to keep the COM of the electrolyte fixed relative to 
the simulation cell. 

Starting from the initial configuration of each simu- 
lated system, an energy minimization was performed, 
followed by a 1 ns NVT equilibration at 300 K for the 
slab systems. For simulations of bulk water ion solvation 
the equilibration process was performed in the NPT en- 
semble, using a weak-coupling pressure regulation with a 
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TABLE I: The system parameters [7 is the H-O-H angle; (j? and Q% x are the permanent molecular dipole and quadrupole 
moments; pi is the liquid density at the center of the slab (estimated error of ±0.005 g/cm 3 )], total interfacial potential <f>i v , th e 



corresponding quadru- {<j>f v ) and dipolar {(j)^,) contributions, and the "isotropic" quadrupolar approximation (<^ cst , eq. 19 1; 
Estimated standard error of ±0.5 mV for the vapor-liquid interface potentials of the studied liquid models (obtained using the 
block averaging method). 



Model 


7(°) 


M °(D) 


QL(DA) 


Pi (g/cm 3 ) 


<f>i v (mV) 








4>? v /4>iv 


SPC 


109.5 


2.347 


8.131 


0.981 


-600.3 


-558.8 


-559.2 


-558.6 


0.932 


SP9 


100.0 


2.347 


5.775 


0.892 


-445.8 


-361.4 


-360.1 


-360.2 


0.808 


SP8 


87.6 


2.347 


3.736 


0.787 


-254.8 


-206.5 


-206.8 


-205.9 


0.811 


SP7 


75.0 


2.347 


2.394 


0.698 


-123.6 


-116.1 


-116.9 


-117.1 


0.946 


SP5 


54.7 


2.347 


1.089 


0.611 


-21.4 


-46.3 


-46.2 


-46.6 




S2N 





2.347 





0.556 


56.8 


0.2 











S2L 





4.065 





0.696 


284.6 


-0.4 











S2D 





2.347 





0.658 


-0.8 


-0.4 












TABLE II: Lennard-Jones parameters a and e, and polarizability a for ions 



Ion 


1(e) 


a (A) 


e (kcal/mol) 


a 


ref. 


Na+ 


1 


2.350 


0.13 


0.24 




F~ 


-1 


3.168 


0.2 


0.974 


|35| 


cr 


-1 


4.339 


0.1 


3.25 


|36| 


Br~ 


-1 


4.700 


0.1 


4.53 


|37| 


r 


-1 


5.150 


0.1 


6.9 


|37| 



target pressure of 1 bar. Subsequently, in both cases, at 
least 5 ns of measurements in the NVT ensemble were 
carried out, using the Berendsen thermostat with the 
configurational degrees of freedom coupled to a heat bath 
with coupling constant t = 1 ps 41j. In the special case 
of polarizability-enabled simulations, the degrees of free- 
dom related to the induced dipole moment of the ion 
were independently coupled to a 1 K heat bath (relax- 
ation time Tdip = 10 ps), ensuring a proper handling of 
the electronic degrees of freedom [42] . All computed pro- 
files spanning the vapor-liquid interface were obtained as 
ensemble averages of the instantaneous profiles evaluated 
in thin slabs (bins) of thickness 0.2 A parallel to the inter- 
face. For the radial quantities measured in the bulk sim- 
ulations of ion solvation, equally distanced, 0.1 A thick, 
spherical shell bins have been employed. Due to the cu- 
bic dimensions of the simulation box the radial profiles 
are relevant up to approximately 15 A. 



V. MD SIMULATION RESULTS 

The MD simulation results show that the bulk region 
density decreases with decreasing molecular quadrupole 
moment (at constant dipole moment), varying by nearly 
a factor of two in going from SPC/E to S2N (table [T). 
The S2N and S2L models form purely dipolar liquids 



with a characteristic chain-like structure arising from the 
head-to-tail alignment of the dipole moments and a lower 
bulk density due to the decrease in hydrogen bond coor- 
dination from four to two. The vapor-water interfacial 
thickness is found to be approximately 3.8 A at 300 K 
for SPC/E and increases along with the slab thickness as 
the molecular quadrupole moment decreases. 



A. Dipolar ordering 

Oricntational (dipolar) ordering of water takes place 
near the 1-v interface, which can be seen by plotting 
P z {z) for the series of model liquids (fig. [2]). Orienta- 
tional double layers were found only for SPC/E and SP9 
with the outer layer dipoles pointing preferentially to- 
wards the vapor phase and in the opposite direction in 
the inner layers (closest to the slab center). For models 
with lower molecular quadrupole moments the molecu- 
lar dipoles point towards the liquid phase. Because of 
the asymmetry created by the oxygen LJ sphere, asym- 
metric purely dipolar liquids (S2N and S2L) still pos- 
sess orientational ordering in the interface region due to 
the hydrophobic forces tending to exclude the oxygen LJ 
sphere from the liquid slab. 
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FIG. 2: (Color on-line) The dipole moment density, P z (z), 
for all studied liquids at the vapor-liquid interface. Dashed 
vertical lines represent the GDS of both interfaces for SPC/E. 



B. Vapor-liquid interface potential: multipole 
contr ibut ions 



In order to illustrate how various multipole moments 
contribute to the vapor-liquid interface potential, we 
obtained the electric potential difference by integration 
of the Poisson equation from the charge density ob- 
tained from the first two terms of the multipole expan- 
sion HEEol, 



d 

dz 



d ( 

dz 



(17) 



leading to the first two (dipolar and quadrupolar) contri- 
butions to the interface potential: 



(z) = 8(j) u {z)+S4^{z) 



dz 



so 



Q zz (z) 

(18) 

since Q zz is taken to vanish in the vapor phase and P z 
vanishes in both bulk (vapor and liquid) phases. For the 
planar interface higher order moments do not contribute. 

The "exact" model interface potential profile and 
the corresponding dipolar and quadrupolar contributions 
(eq. 18 1 obtained directly from the simulation data (us- 
ing, respectively, the "exact" partial charge density, p, 
and the multipole contributions of eq. 17) are illus- 
trated in fig. [3j The total vapor-liquid potential reaches 
— 600mV for SPC/E, in agreement with previous val- 
ues [5] , but decreases in amplitude with decreasing molec- 
ular quadrupole moment (SP9 to SP5). The asymmetric 
purely dipolar liquids, on the other hand, have positive 
interface potentials, whereas, as expected by symmetry 
considerations, the fully symmetric model S2D gives a 
null result (table [I]) . The two components of the inter- 
face potential reveal very different types of profiles with 
the quadrupolar potential being negative, as expected 



n -0.2 




FIG. 3: (Color on-line) (a) Total, (b) dipolar and (c) 
quadrupolar potential profiles, 8<j){z) (in volts) for the stud- 
ied liquids characterized by positive molecular quadrupolar 
moments. 



for models with positive molecular quadrupole moments. 
For the SPC/E model the quadrupolar contribution rep- 
resents more than 90% of the total. This contribution 
decreases rapidly with decreasing molecular quadrupo- 
lar moment and becomes comparable in absolute value 
with the (positive) dipolar contribution for SP5, the near 
cancelation in this case leading to a very low (negative) 
total value. For SPC/E, SP9-SP7 the quadrupole contri- 
bution provides the major contribution to the interface 
potential (figs. [3] and |4j table [B and thus to the large 
interfacial electric fields (~ 1 V/nm for SPC/E) directed 
towards the liquid phase over a substantial part of the 
interfacial region For SPC/E and SP9 the inner- 

most dipoles tend to follow this field, creating in turn an 
opposing dipolar field that acts to align the outermost 
molecular dipoles in the opposite direction. Although 
the quadrupolar potential profile is always monotonic, 
the dipolar one shows a minimum close to the Gibbs di- 
viding surface (GDS) for both SPC/E and SP9 models 
due to the dipolar orientational bilayers. These results re- 
veal a subtle interplay between the dominant quadrupo- 
lar contribution and the dipolar response. 
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FIG. 4: (Color on-line) Vapor-liquid interface potential drops 
for all studied models: standard (blue) and nipped (F)-charge 
(black) with respect to \ Q% X \ (molecular quadrupole moment). 



Even if the system is not isotropic in the interfacial 
region, it is possible to generalize the approach pre- 
sented in [5J [TU1 H3] to construct a simple but extremely 
accurate "isotropic" approximation for the quadrupolar 
contribution to the local vapor-liquid interface potential 
using only the water density profile and the molecular 
quadrupole moment Q° evaluated in a local reference 
frame with the y-axis along the dipole vector and the 
z-axis out of the molecular plane: 



b ^ z >- 6£n 3 ' 



(19) 



where c(z) is the local liquid number density taken from 
the simulations. In this reference frame the only non-zero 
component is Q xx and therefore in this case Tr Q° = Q xx . 
The estimate (f>f v cst — 5<fif st (zi) for the quadrupolar con- 
tribution is in excellent agreement with direct determi- 
nations (table [l| . We have also checked that this sim- 
ple quadrupolar estimation [19] provides a very good ap- 
proximation to the full oscillatory membrane-water sur- 
face potential |21) . confirming that the air- water interfa- 
cial and membrane-water surface potentials are mainly 
determined by the local water density and molecular 
quadrupole moment. Furthermore, we propose that this 
method can be used to estimate the quadrupolar poten- 
tial contribution of any liquid, irrespective of its bulk 
molecular quadrupole moment, even those obtained from 
ab initio quantum mechanical calculations of liquid wa- 
ter. As an illustration, we have checked that when ab ini- 
tio values for Tr Q° [5?1|35], are injected into the simple 



approximation Eq. 19 we find quadrupolar potentials in 
reasonable agreement with the (quadrupole dominated) 
total interface potentials (+3 eV) extracted directly from 
the quantum mechanical calculations [T2J [T3] . 



C. Solute-Liquid interface potential 

In the presence of a solute the planar vapor-liquid in- 
terface potential has as counterpart the microscopic po- 
tential between the solute center and the surrounding 
liquid (of which the first two multipole terms can be ob- 
tained from the microscopic analog of eq. 18 29J ) . 



The local radial solute-liquid potential profile, defined 



as 



Scj) r (r) = 4> r (r) - 4> r (r s ) 



(20) 



[with 4> r (r s ) = at the center of the solute r s = 0] , is ob- 
tained from the radial charge density p r (r) by integrating 
the Poisson equation pi in spherical coordinates: 



b r {r) =- J E r (r')dr' = 

t> 



Pr (r') (r'f ( 1 1 



so 



dr'. 



(21) 

The dipolar component, 8<f>® (r), is determined from the 
dipole moment density P r (r) as: 



6?(r) 



dr' 



(22) 



The corresponding radial electric fields, E r (r) and 
E® (r), are determined from equations (10 1 and (15). 



The quadrupolar contribution to the total radial in- 
terface potential, 5<fiQ (r) is accessible from the simula- 
tion data via the radial dependence of the quadrupole 
moment density written in spherical coordinates, Q (r). 
We begin by writing the quadrupolar Poisson equation 
(13) in Cartesian coordinates (centered at the solute po- 
sition r s — 0), with the tensor elements of the Cartesian 
quadrupole moment density, Q a p: 



(23) 



ct,/3 



Using the divergence theorem, we find: 



£ $E Q -dS= ©X^^-r^dS. (24) 

s s a >P 

Letting F 'g = ^ ^ a Qaj3 and S — Airr 2 , we obtain the ra- 

dial dependence of the quadrupolar electric field E® (r), 
by integrating over the angular degrees of freedom: 



AttsoEQ (r) = © F r sin 

5' 



(25) 



with F r = e r -Ffj. After performing the angular integrals, 
we obtain the radial dependence of the quadrupolar elec- 
tric field: 

e E? (r) = ^ + - (3Q' rr - TrQ) . (26) 
or r 
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Since (r) = -^§^ and ^ (0) = 0, we arrive at 
the final solution for the local quadrupolar potential at 
position r, with respect to the center of the solute, 
(r) — </>5? (r s ), in terms of two components: 



with 



Stf' 1 (r) 
Stf> 2 (r) 



Q'rr (r) 



(27) 



(28) 



dr 



£o J r 
o 



r [TrQ(r')-3Q; r (r')]. (29) 



The second contribution is generated by the symmetry 
breaking of the diagonal components of Q' in the solute- 
liquid interfacial region. The detailed calculations lead- 
ing from Q (r) to 5(j)Q (r) will be presented elsewhere [53] . 
Far from the solute center the radial solute-liquid inter- 
face potential and the various multipole components tend 
to their respective asymptotic values: <j)i s = Scf> r (oo), 
0% = 5<tf?(oo), and cf>? s = S^(oo). 

Because the solute-liquid interface is curved, the as- 
sociated potential depends on the size of the cavity and 
is therefore ion specific. Thus, there is a first potential 
drop when going from the vapor into the liquid phase, 
followed by an overall increase near the solute, yielding 
a smaller, overall negative, vapor-liquid-solute potential 
drop, (j) sv (fig. [lj. To obtain a fuller picture, we have 
extracted the solute-liquid interfacial potentials (and the 
dipolar/quadrupolar components) from MD simulations 
of neutral ion-like solutes, fixed at the center of a cu- 



bic box of bulk SPC/E water (using eqs. |2lj [22| and[27|. 
The solute-liquid potential 4>i s and the corresponding 
dipolar contribution ipf^ are obtained, respectively, from 
the radial charge density and the dipole moment density 
distributions. The higher order multipole contributions, 
(f>i s — 4>f s , are dominant and, due to Is interface curvature 
effects, not only is the amplitude of the Is quadrupo- 
lar contribution different from the Iv one, but multipole 
terms beyond the quadrupolar one play a non-negligible 
role (tables |III| and IV I [29 . For the halide-like neutral 
solutes 4>i s decreases in amplitude with increasing solute 
size and is smaller than for the neutral sodium like solute, 
Na°. For 1°, 4>i s increases in amplitude by about 6% when 
the polarizability is turned on. Although the (ft®' 1 con- 
tribution is dominant in (f>i s , it simply serves, as we shall 
see below, to cancel the large vapor-liquid quadrupolar 
contribution to the direct contribution to the free energy 
of ion solvation, AGo- The dipole contribution (f)f s goes 
from negative values for small halide-like solutes (F°) to 
positive values for larger ones. 



AGo = q(4>iv — 4>is) (eq. [3]) , dominated by the difference 
between the quadrupolar Iv and Is contributions, which 
do not cancel due to strong curvature effects for the small 
ions under study. Because the cpf' 1 component depends 
only on the bulk solvent properties, it is in principle equal 
to 4>® v and therefore the two terms should cancel in AGo, 
leaving the domi nant ion spe cific quadrupolar contribu- 



tion, 



(tables 



III 



and IVb . The dipolar contribution 
and multipolar contributions higher than quadrupolar 
play non-negligible, but secondary roles, in AGo- Our 
results for polarizable I~ also reveal that ionic polar- 
izability plays a minor but non-negligible role in deter- 
mining AGq (tables III and IV). The key point is that 



AGo is much smaller in amplitude than the simple direct 
estimate (|AG | = 0.600 eV) because of the strong par- 
tial cancelation of the two interface potentials, <f>i v and 
(pi s . Our results show that the direct contribution, AGo 
is reduced to 30-40% of the simple direct estimate AG 
and depends on the sign of the ion charge in such a way 
as to favor the solvation of cations Na + with respect to 
the anions (halogens). For the halogens AGo is positive 
and increases with increasing ion size, thus augmenting 
the interface propensity of large anions like I - . Because 
the amplitude of AGo (~ 0.2 eV) is still between 5 and 
10% of the dominant polarization contribution, AG po i, 
it is on par with other important contributions (such as 
the hydrophobic one [23 123 123 123 E5]) and therefore 
must be taken into account correctly when considering 
ion-specific effects. 

Because of charge-dipole and charge-quadrupole cou- 
pling, when the ion charge is "turned on" the solvent 
molecules reorient themselves in order to accommodate 
the solute, optimizing as much as possible the orienta- 
tion and number of hydrogen bonds: the ion polarizes 
the medium around it and induces a radial potential 
difference, 4>\° s n , dominated by the "long range" dipo- 



lar contribution. The difference between (/>J° n and <pi s 
extracted from the simulations then determines AG po i 
(eq.[4j. With our choice of effective ion radius, the simple 
mesoscopic Born approximation AGp ol is in reasonably 
good agreement with the microscopic polarization contri- 
bution, AGpoi (table [V| . The polarization contribution 
to the solvation free energy, which is the main effect fa- 
voring high ion solvation, decreases in amplitude with 
increasing halogen size. We note also that the simula- 
tion results for the electrostatic contribution to the ionic 
free energy of solvation, AG^g, follow the experimen- 
tal trend, AG™, reasonably well (table [vj), despite the 
neglect of certain entropic and enthalpic contributions 
(arising, for example, from the short range repulsion and 
the long range van der Waals attraction). 



D. Free energy of ion solvation 

The solute-liquid potential was then used to evaluate 
the direct contribution to the free energy of ion solvation, 



VI. CONCLUSION 

We have studied a series of liquid models that interpo- 
late between SPC/E water and pure dipolar liquids and 
shown that the quadrupolar component of the vapor- 
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TABLE III: Solute-liquid interface potentials with multipolar contributions for neutral solutes solvated in SPC/E water (for 
which the vapor-liquid interface potential is <f>i u = —600.3 mV with 4>f v = —559.2 mV and <j>% = —41.5 mV). In principle 
^is 1 = 4>i v - The less than 2% difference between the values obtained from the water slab and bulk simulations can be 
attributed to the differences in densities generated by the barostat used in the NPT bulk simulations (cf. eq 19 1. To directly 
compare solute-liquid and vapor-liquid interface potentials, we correct for this small systematic liquid density difference by 
defining a solute-dependent rescaled vapor-liquid interface potential, (f>iv — C<f>i v , where C = 0^'V<^ (rescaled multipole 
vapor-liquid interface components are defined similarly). The estimated standard error for the interface potentials is 0.15 mV 
and pol denotes polarizable. 



Solute 


4W<^ (%) 


<MmV) 


4>u (mV) 


<t>l (mV) 


tf,' 1 K) 


</>?/ (mV) ) 




Na° 


66.7 


-405.7 


-78.8 


-423.6 


-566.7 


143.1 


96.7 


F° 


62.4 


-380.0 


-41.6 


-415.3 


-567.3 


152.0 


76.9 


Cl° 


60.3 


-369.0 


0.3 


-426.8 


-569.9 


143.1 


57.5 


Br° 


58.8 


-359.8 


17.8 


-430.8 


-569.6 


138.8 


53.2 


1° 


58.4 


-357.5 


32.3 


-437.6 


-570.5 


132.9 


47.8 


1° (pol) 


62.1 


-380.1 


14.4 


-440.0 


-570.5 


130.5 


45.5 



TABLE IV: Comparison of the different multipole contributions to the direct electrostatic solvation-free energies, AGS = 
q(4>* v — <t>is), obtained using the rescaled vapor-liquid interface potentials (*), of various ions as obtained directly from the 
simulations (see table [n} (AG' ~ ±0.6003 eV); estimated errors: 0.002 eV 



Solute 


AGS (eV) 


Dipole* (%) 


Quadrupole* (%) 


Higher order* (%) 


Na+ 


-0.2030 


-18.1 


70.5 


47.6 


F" 


0.2294 


0.2 


66.3 


33.5 


cr 


0.2432 


17.5 


58.8 


23.6 


Br~ 


0.2521 


23.8 


55.1 


21.1 


I" 


0.2554 


29.2 


52.0 


18.8 


I- (pol) 


0.2327 


24.4 


56.0 


19.6 



TABLE V: Comparison of electrostatic solvation-free energies for various ions as obtained directly from the simulations, AGgg = 
AGo + AGpoi, with the estimate AG^est — AG ± AGp 0l obtained from the truncated direct contribution [AG ~ ±0.6 eV 
for monovalent anions (resp. cations)] and the Born approximation, AGp 0l . Ri is chosen as the distance from the ion center 
where the ion- water radial distribution functions first reach 1 [29] . To correct for the systematic differences in liquid density 
between water slab and bulk simulations, the rescaled vapor-liquid interface potentials, </>*„, are used in the evaluation of the 
direct contribution: AGS = l(4>iv ~ 4>is) (see table Inl. Estimated error of 0.025A for Ri; Estimated standard errors: 0.002 eV 
for AGo; 0.05 eV for the other solvation-free energies. (1 eV= 23.06 kcal/mol~ 40fcBT) 



Solute 


R l (A) 


AGS (eV) 


AG^, (eV) 


AG P oi (eV) 


AGf|, e8t (eV) 


AG [ i§ (eV) 


AGi° x n P (eV) gg 


Na" 


1.05 


0.2030 


-6.76 


-6.88 


-6.57 


-6.68 




F~ 


1.55 


0.2294 


-4.58 


-4.46 


-4.36 


-4.23 


-4.50 


cr 


2.00 


0.2432 


-3.55 


-3.34 


-3.32 


-3.10 


-3.64 


Br" 


2.25 


0.2521 


-3.15 


-3.01 


-2.91 


-2.76 


-3.30 


r 


2.45 


0.2554 


-2.90 


-2.63 


-2.66 


-2.38 


-2.90 


Na+ 


2.20 


-0.2030 


-3.23 


-3.88 


-3.42 


-4.08 


-4.26 


1+ 


3.55 


-0.2554 


-2.00 


-1.79 


-2.24 


-2.04 





liquid interfacial potential typically dominates for the 
studied liquids possessing a non-zero quadrupolar mo- 
ment. In an effort to elucidate the different ion-specific 
contributions to the free energy of solvation, we have 



shed light on the key role played by the solute-liquid 
interface potential and demonstrated that it leads to a 
strong reduction in the direct electrostatic contribution 
with respect to previous estimates based solely on the 
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vapor-liquid potential. 

We propose that the same mechanism would be at 
play if the point partial charge distribution of the sol- 
vent extracted from classical MD simulations were re- 
placed by the more realistic extended charge distribu- 
tions found in ab initio calculations. Indeed, a coarse 
graining procedure for the electric potential proposed re- 
cently [13] . which corrects for regions inaccessible to ionic 
probes, shows that, encouragingly, both ab initio and 
point charge coarse grained potentials converge to val- 
ues that are compatible with the results for AG = q<f> S v 
presented in table |Vj Finally the dominant electrostatic 
polarization contribution to the free energy of solvation 
was found to agree reasonably well with a Born-type ap- 
proximation. We conclude that the direct interface po- 
tential contribution to the ionic free energy of solvation 
(or PMF) can neither be estimated using the point ion 
approximation (leading to a gross overestimate), nor be 
neglected entirely - the two approximations commonly 
adopted in the current literature. One important corol- 
lary that can be drawn from our study is that even purely 
quadrupolar liquids should give rise to an interface po- 
tential contribution to the ionic solvation free energy be- 
cause of the incomplete cancelation of the to and Is com- 
ponents (due to solute curvature effects). Another impor- 
tant corollary is that dipolar statistical mechanical mod- 
els of water near interfaces [47) , despite their evident in- 
trinsic interest, may be inadequate for capturing certain 
effects arising from interfacial potentials, because of the 
seemingly essential role played by the water quadrupole 
moment. 

The simulation results presented here should be not 
only a useful building block in the quest for construct- 
ing the physically relevant mesoscopic components of the 
ion free energies of solvation (or PMF) , but also provide 
a challenge for statistical mechanical approaches used to 
analyze the subtle interplay between short range steric, 
dipolar, and quadrupolar interactions in determining the 
properties of polar fluids (in particular the dipole distri- 
bution near interfaces) and their influence on ions. 



An important outcome of a reliable mesoscopic ap- 
proach would be a greatly enhanced comprehension of 
the underlying physics of ion distributions in inhomoge- 
neous dielectric settings with important applications in 
colloidal science, nanotechnology (ion transport in arti- 
ficial nano-pores, or nano- filtration) , and biophysics (ion 
channels, biological membranes, DNA) [15-18]. We also 
expect that the results presented here transcend the par- 
ticular chosen models and thus qualitatively illustrate im- 
portant physicochemical mechanisms at play in ion par- 
titioning within inhomogeneous dielectric media. 

We are currently extending the approach developed 
here for the ionic solvation free-energy to the study of 
the local ionic PMF near aqueous interfaces and surfaces 
with and without the potentially important effect of wa- 
ter and ion polarizability [2"TH2"3l 15211431 |4"5] . 
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